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ABSTRACT 



This research designed and simulated an adaptive attitude control system for 
the Crew Equipment/Retriever ( CER ) during autonomous attitude hold and large 
angle or slewing maneuvers. The CER is a proposed space robot that deploys 
from the Space Station and retrieves any lost equipment or incapacitated 
astronauts. The moment of inertia tensor for the CER and acquired target is not 
known a priori. In this research, the moment of inertia tensor is estimated by a 
Kalman filter and used to update the derived linear quadratic regulator ( LQR ) 
and quaternion feedback regulator ( QFR ) control laws. Computer simulation 
results show that during attitude hold the adaptive LQR design stabilizes the CER 
and provides a more fuel efficient controller effort: as compared with a 

previously designed nonadaptive minimum time controller and a nonadaptive 
LQR design. Computer simulation results of slewing maneuvers show that the 
adaptive QFR design provides a more fuel efficient controller: as compared with 
a nonadaptive QFR design. 
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I. INTRODUCTION 



A. CERS CONCEPT DESCRIPTION 
1. CERS Origin and Purpose 

The National Aeronautics and Space Administration (NASA), Johnson 
Space Center Space Station Projects Office sent out a Request for Proposal in 
May 1987 as part of Space Station Work Package 2. This Request for Proposal , 
including an added Amendment 7, defined a requirement to provide for the 
capability to rescue an incapacitated external-vehicular activity crewman and to 
retrieve equipment inadvertently detached from the Space Station. [Ref.l: p. L- 
2- 14a] 

McDonnell Douglas Astronautics Company (MDAC) responded to this 
Request for Proposal in September 1987 with a practical, low cost retriever 
concept. This concept was referred to as the Crew/Equipment Retrieval System 
(CERS). [Ref. 2: p. 1] This overall system consisted of a crew and equipment 
retriever vehicle (CER) and other Space Station based support systems. 

The overall mission of the CER is to: deploy from the Space Station, 
acquire and capture the designated target, and return to the Space Station. A 
summary of the CERS capabilities, as defined by MDAC, is listed as follows: 
[Ref. 2: p. 9] 

1. Retrieve an 850 pound target (includes 10% safety margin); 

2. Total Deployment time of 120 minutes; 

3. Retriever activated and deployed without assistance from 
an external-vehicular activity crewman; 

4. Retriever senses own attitude, range, and range to target; 
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5. Retriever can be remotely operated from the Space Station; 

6. Retriever accommodates a worst-case separation speed of 3.5 ft/sec; 

7. Retriever senses and controls its own attitude with and without a target; 

8. Retriever has attitude hold and three-axis translation capability. 

2. CERS Baseline Configuration 

Figure 1 shows the CER and its Space Station support systems. More 
detailed descriptions of proposed hardware and software are contained in Ref. 2 

(pp. 20-68). 

A simple representation of the CER for attitude dynamics analysis was 
developed in Ref. 3 and is shown in Figure 2. The characteristics of the baseline 
configuration developed by MDAC are listed as follows: [Ref. 2: p. 24] 

1 . 850 pounds total weight; 

2. Three-axis (six degrees of freedom) stabilized; 

3. Remote tele-operated free flyer; 

4. Use of 24 cold Nitrogen (N 2 ) jet thrusters rated at 1.0 lbf each; 

5. Attitude control is accomplished by firing thrusters in pairs; 

6. Maximum control torques: 

Roll Axis 3 ft-lbf; 

Pitch Axis 3 ft-lbf; 

Yaw Axis 4 ft-lbf. 
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Figure 1. CERS Major Components [Ref. 2: p. 15] 
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Figure 2. CER Baseline Configuration [Ref. 3: p. 8] 
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An important feature of the CER is its ability to control its attitude. After 
acquiring a target the CER must stabilize its attitude and perform an attitude 
reorientation prior to returning to the Space Station. 

The CER’s attitude control problem is very different from the norm. 
Most spacecraft are designed such that any changes in their moments of inertia 
are minimized. The control devices are, moreover, placed so as to act along the 
principal axes of the body. This positioning minimizes any gyroscopic coupling 
between torque applied to any one axis and rotation about another axis. Neither 
of these two conditions hold true for the CER after it acquires a target since the 
target can be as massive as the CER by itself. 

B. THESIS OBJECTIVES 

The overall objective of this thesis is to design adaptive attitude control laws 
for the CER with and without a target during large angle or slewing maneuvers 
and during autonomous attitude hold for all mission phases. Previous thesis 
research [Ref. 3] modeled the CER and designed time-optimal and weighted 
time-fuel optimal single-axis control systems and tested these control systems by 
applying them to several worst case target scenarios. This research analyzes the 
complete nonlinear three-axis control problem for small angle motion or attitude 
hold and large angle or slewing motion. Attitude hold pertains to attitude 
control in the presence of small disturbances while slewing motion pertains to 
attitude reorientation. The control law designs are adaptive in that a key system 
parameter, the moment of inertia of the CER and acquired target, is not known a 
priori and must be estimated. A subsidiary goal of the slewing motion control 
law is to accomplish this reorientation in an optimal fashion: an eigenaxis 
rotation. Both control laws must be able to deal with a non-diagonal moment of 
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inertia tensor since after target acquisition, the thrusters no longer act along the 
principal axes of the body. 

C. THESIS ORGANIZATION 

In Chapter II, general spacecraft attitude kinematics and dynamics are 
developed. These equations of motion are then applied to the CER. The moment 
of inertia tensors for the CER with and without a target are calculated. 

Chapter III derives two control law designs. One design is developed for 
attitude hold while another design is developed for large angle or slewing 
maneuvers. Central to each control law design is the knowledge of the CER 
moment of inertia tensor. 

In Chapter IV, an estimation scheme is developed that provides the above 
mentioned control laws with an estimate of the CER moment of inertia tensor. 
This estimation scheme is based on a rather unusual application of the Kalman 
Filter. 

Chapter V presents the computer simulation results of applying each control 
law and estimation scheme. Control system design issues and practical 
implementation details are also discussed. 

Conclusions based on the computer simulation results are presented in 
Chapter VI. In addition, recommendations for future research are discussed. 
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II. ATTITUDE KINEMATICS AND DYNAMICS 



The equations of motion for any rotating rigid body can be divided into two 
sets: the Kinematic Equations of Motion and the Dynamic Equations of Motion. 
Kinematics studies motion without considering the forces that cause that motion. 
The Kinematic Equations of Motion are a set of first order differential equations 
that specify the time evolution of the chosen attitude parameters. The Dynamic 
Equations of Motion , meanwhile, take into account the forces that cause 
rotational motion and express the time evolution of the angular velocity of the 
rigid body. [Ref. 4: p. 510] 

A. ROTATIONAL KINEMATICS 

1. Direction Cosine Matrix 

Any general vector r can be written in terms of its magnitude and 
direction. The direction can be represented as a unit vector referenced to some 
previously defined coordinate or reference frame as shown in Figure 3. This 
unit vector is made up of components known as direction cosines : [Ref. 5: p. 9] 

r = rr = r[(cosa)n 1 +(cosP)n 2 + (cosy)h 3 ] (2.1) 

where r is a scalar magnitude, and r is a unit vector whose components are 
referenced to the three orthogonal axes of the reference frame ii. Note the 
following symbolic conventions used throughout this thesis: 

1. The underline bar denotes a vector, r 

2. The hat symbol denotes a unit vector, r; 

3. The double underline bar denotes a matrix, W. 
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ii =rcosa 
r 2 = rcosP 
r 3 = rcosy 


r = if = r[(cosa)n! + (cosp)n 2 + (cosy)n 3 ] 



Figure 3. Direction Cosines [Ref. 5: p. 9] 
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Now consider a unit vector b with components along the orthogonal 
body axes of a rigid body, i.e., a spacecraft as shown in Figure 4. Let n be a 
unit vector with components along three orthogonal directions that are fixed in 
space. These two unit vectors are related by the following transformation: [Ref. 
5: p. 9] 



b = Tn (2.2) 

The transformation is defined by T, the 3x3 direction cosines matrix (DCM). 
This DCM is critical to the field of spacecraft dynamics and control and Ref. 5 
addresses several of its important properties. The most important of these 
properties are summarized below: 

1. A DCM exists for any pair of orthogonal sets of three axes; 

2. The DCM is an orthogonal matrix and its inverse equals its transpose; 

3. A DCM can be built up from successive rotations about the axes. 

The last property is best explained by an example. Define a sequence of 
reference frames related by the following transformations: 



II 

':=> 


(2.3) 


b = T 2 a ; 


(2.4) 


c — T 3 b, 


(2.5) 



and therefore the reference frames c and n are related by: 
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Figure 4. Spacecraft Body Axes 
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( 2 . 6 ) 



c = T 3 T 2 T 1 ri. 

In other words, one overall DCM can be formed as a product of DCMs: 

t 321 =t 3 t 2 t;. (2.7) 



2. Euler Angles 

The relative orientation of two orthogonal reference frames can be 
defined in terms of three angles [Ref. 5: p. 16]. This idea was first introduced 
by Euler in the early eighteenth century and is synonymous with the idea of 
parameterizing the previously discussed nine element DCM with only three 
independent parameters. 

The classical Euler angles, for which there are twelve distinct cases, 
define an arbitrary orientation by using the successive rotational transformation 
property of the DCM. That is, a sequence of three elementary rigid right handed 
rotations about instantaneously fixed axes is used to build up an overall DCM that 
represents the transformation from one orientation to a different orientation. 
[Ref. 5: p. 17] 

This thesis employs the 3-2-1 or yaw -pitch- roll Euler angle sequence 
since it is most commonly used in aircraft and spacecraft applications. This 
sequence is produced by initially lining up both the body axes and the fixed or 
inertial axes. A rotation about the z or number three axis is then performed, 
producing a new y and x axis . A second rotation about the new y or number 
two axis is then enacted. Finally, a rotation about the new x axis finishes the 
rotation sequence. Figure 5 depicts the formulation of this sequence. 
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Body Axes and Inertial Axes 
are initially lined up 



A rotation is then performed 





Next, a rotation about the 
y or pitch axis is performed 



Finally, a rotation about the 
x or roll axis is performed 





The overall transformation is 



T321 ~ TiT 2 T 3 
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Figure 5. The 3-2-1 Euler Angles 
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Reference 5 illustrates the general process for obtaining the kinematic 
differential equations for a chosen set of Euler angles. Wertz [Ref. 4: p. 765] 
lists the kinematic equations of motion for the twelve possible Euler angle 
representations. The kinematic differential equations for the chosen 3-2-1 Euler 
angle set are as follows: 



\j f = (co Y sin(<p) + co z cos(<p)) / cos(0) ; 


(2.8) 


0 = co Y cos(<p) - co z sin(cp) ; 


(2.9) 


(p = co x + (co Y sin(cp) + co z cos(cp)) tan(0) ; 


(2.10) 



where: 

1. \|/ is the yaw angle; 

2. 0 is the pitch angle; 

3. cp is the roll angle; 

4. The vector co is the angular velocity vector and is composed of components 
along each of the body axes. 

3. Euler’s Principal Rotation Theorem 

Junkins [Ref. 5: p. 26] states that Euler is generally credited with being 

responsible for the Principal Rotation Theorem : 



A rigid body can be brought from an arbitrary initial 
orientation to an arbitrary final orientation by a single rotation of 
the body through a principle angle about a principal line; the 
principal line being a judicious axis fixed in the body and fixed in 
space. 
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This concept, displayed in Figure 6, allows the DCM to be 
parameterized in terms of the principal angle (|> and principal line £. In 
mathematical terms, the principal line corresponds to the eigenvector of the 
DCM: for the eigenvalue + 1. Therefore, given any DCM one can solve for the 
principal line and angle and reduce the general angular displacement to a single 
rotation about a fixed line. [Ref. 5: p. 27] 

4. Euler Parameters 

In conjunction with the Principal Rotation Theorem , Euler defined four 
parameters in terms of the principal line and principal angle. These Euler 
Parameters are as follows: 



Po = cos(<J) / 2) (2.11) 

Pi = £ x sin(<{) / 2) 

P 2 = £ 2 sin(<|> / 2) 

P 3 =^ 3 sin«J>/2) 

where t l9 £ 2 , £ 3 are the components of the unit vector along the principal line £ 
and (j) is the principal angle. 

The DCM is therefore parameterized in terms of the above Euler 
Parameters and this allows the relative orientation of two orthogonal reference 
frames to be represented by four parameters. One of these parameters is 
redundant since the DCM was previously shown to be parameterized by three 
Euler angles. This redundancy manifests itself in the following constraint 
equation: 
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Figure 6. Euler's Principal Rotation Theorem [Ref. 5: p. 25] 
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Pc 2 +P, 2 +P 2 2 +P3 2 =l- 



( 2 . 12 ) 



Various algorithms have been developed that determine the Euler 
parameters from a given DCM and convert back and forth between Euler angles 
and parameters. These algorithms have been used extensively in the computer 
simulation programs for this thesis and are included in the Appendix. [Ref. 5: 
pp. 31-35] 

By differentiating the inverse relationships between the Euler 
parameters and the elements of the DCM, and making a few substitutions, the 
kinematic differential equations in terms of the Euler parameters can be 
formulated as: [Ref. 5: p. 35] 



'Po' 




'0 


-CO, 


-co 2 


-co 3 " 


'Pc' 


p, 


_ J_ 


CDi 


0 


co 3 


-co 2 


Pi 


p 2 


"2 


co 2 


-co 3 


0 


COi 


P 2 


A. 




_co 3 


co 2 


-COi 


0 


,p3. 



(2.13) 



Note that throughout the literature on attitude dynamics and control, the 
Euler parameters are sometimes referred to as quaternions and are formulated as 
follows: [Ref. 4: p. 414] 



q, = /j sin((J) / 2) (2.14) 

q 2 = l 2 sin(<j> / 2) 
q 3 = £ 3 sin((J) / 2) 
q 4 = cos((J) / 2) 
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The constraint equation (2.12) is also applicable and a conversion from one 
representation of Euler parameters to this quaternion representation allows the 
following identification: 



q 4 s Po (2.15) 

q, - Pi> i = 1» 2, 3. 

As seen above in equation (2.15), the difference between these two 
representations is rather insignificant. Note that this author has chosen to use the 
Euler parameter representation of equation (2.11) but throughout this thesis the 
terms Euler parameters and quaternions are used interchangeably. 

5. Parameterization Discussion 

The previous sections briefly demonstrate that there exist several 
choices when representing the orientation of a rigid body. Euler angles are 
easier to visualize and are more popular but suffer from one large draw back: 
the presence of mathematical singularities at certain angles. Equation (2.8), for 
example, experiences a singularity when the cosine of the pitch angle goes to 
zero. This can be a real problem in terms of numerical computations during 
computer simulation or software running control system algorithms. The use of 
Euler parameters eliminates the use of trigonometric functions and their 
singularities but they are harder to visualize. 

Therefore, the intended application should dictate the choice of 
kinematic differential equations. In this thesis, the 3-2-1 Euler angle set is used 
to define the kinematics for small angle motion during attitude hold. For large 
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angle or slewing motion and control, the less widely used but much more 
practical Euler parameters are used to represent the kinematics. 

B. ROTATIONAL DYNAMICS 

1. General Rigid Body 

The Eulerian Rotational Equations of Motion for a rigid body subject to 
applied control torques are well known and are typically represented as: [Ref. 5: 
pp. 49-52] 



cb = -I” 1 coIcd + 1” 1 u 



(2.16) 



where: 

1 . co is the angular acceleration vector; 

2. I is the moment of inertia tensor for the rigid body; 

3. co is the angular velocity vector; 

4. u is the vector of applied control torques; 

5. cd is the skew symmetric angular velocity matrix defined as: 



0 


— cd 3 


cd 2 


cd 3 


0 


-CD! 


_cd 2 


CD t 


0 



(2.17) 



2. The CER 

The dynamics of the CER can be represented by equation (2.16) 
provided that the rigid body assumption is used. The only system parameter in 
equation (2.16) that is specific to the CER is its moment of inertia tensor. This 
moment of inertia tensor will change whenever the CER captures a target. 
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Hansen [Ref. 3] calculated the moment of inertia for the CER without 
target by assuming an 850 pound total system weight symmetrically distributed 
about the center of gravity. The resulting matrix is diagonal in the CER's body 
coordinate system and has the units of slug-foot 2 : 



"39.6 

0 

0 



0 

55 

0 



0 ‘ 

0 

55 



(2.18) 



The total moment of inertia of the CER with a target, as previously 
discussed, is not known a priori since it depends on the specific mass and 
moment of inertia of the acquired target. Hansen [Ref. 3: p. 17] utilized the 
parallel axis theorem to calculate the total moment of inertia for the CER and an 
acquired target 



M t M 2 2 
M, + M 2 — 



(2.19) 



where Mj is the mass of the CER, M 2 is the mass of the target, and R 2 is the 
skew symmetric matrix derived from the vector r 2 between the center of gravity 
of the CER and the target: 



R,= 



0 

r z 




( 2 . 20 ) 



The fact that vectors can only be added together if they are expressed in the same 
coordinate system also applies to the moments of inertia in equation (2.19). Here 
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again, the DCM T plays the role of transforming a moment of inertia from one 
reference frame to another reference frame as follows: 



l TARGET 




( 2 . 21 ) 



where the superscript T refers to taking the transpose of the given matrix, the 
superscripts b and n refer to the reference frame, and I is the moment of inertia 

tensor. 

Hansen [Ref. 3: pp. 17-21] further utilized equations (2.18-2.21) to 
calculate five worst case target scenarios, each with an associated total moment of 
inertia tensor. This thesis used the Case Two scenario most frequently for 
control system design analysis and computer simulation. This scenario 
corresponds to an astronaut with manned maneuvering unit in the target net and 
has a total moment of inertia ( slug-foot 2 ) of: 



C. STATE VARIABLE REPRESENTATION 

The total set of equations representing the CER attitude kinematics and 
dynamics is composed of equation (2.16) and either equations (2.8-2.10) for the 
3-2-1 Euler angle representation or equation (2.13) for the Euler parameter 
representation. These equations are nonlinear in nature but for initial control 
system design and as an approximation when small angle motion is considered, 
equation (2.16) and equations (2.8-2.10) can be approximated as: 



L 



1 TOTAL 



112.9 2.4 -111.9' 

2.4 534.9 6.4 

-111.9 6.4 497.6 



( 2 . 22 ) 
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where X = [v 0 <p co z co y co x ] T , u is the applied control torques and: 





'0 


0 


0‘ 




0 


0 


0 


B = 


0 


0 


0 




In 


Itt 


1 13 




I* 


la 


£ 




In 


I* 







'0 


0 


0 


1 


0 


O' 




0 


0 


0 


0 


1 


0 


A = 


0 


0 


0 


0 


0 


1 




0 


0 


0 


0 


0 


0 




0 


0 


0 


0 


0 


0 




0 


0 


0 


0 


0 


0 



(2.24) 



(2.25) 



In order to complete the state space equations, the plant defined by equations 
(2.23-2.25) must be accompanied by an output equation: 

Y = CX + Du. (2.26) 

This thesis assumes that all the states are measurable and there is no direct 
coupling between control input and the state outputs. Therefore, the following 
definitions are made: 



"1 0 0 0 0 0 " 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0 ’ 

0 0 0 0 1 0 

0 0 0 0 0 1 



(2.27) 



"0 0 0 " 
0 0 0 
0 0 0 
0 0 0 
0 0 0 
0 0 0 



(2.28) 



Note that the C matrix is simply the identity matrix and the D matrix is simply a 
matrix of zeros. Therefore, equation (2.26) can be reduced to: 

Y = X. (2.27) 

Realistically, equation (2.27) should contain some added errors or noise due to 
the fact that sensors have limited accuracy and do introduce noise into any actual 
system. A more detailed discussion on sensors is contained in Chapter III and a 
related discussion on computer simulation implementation details is found in 
Chapter V. 
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III. ATTITUDE CONTROL LAW DESIGN 



The development of the attitude equations of motion in the previous chapter 
and the design of control laws or algorithms in this chapter are based on one 
fundamental assumption: attitude motion of a spacecraft can be approximately 
decoupled from its orbital motion. For the purpose of attitude control design, 
therefore, the spacecraft is almost universally considered to have only rotational 
degrees of freedom about its center of mass which is fixed to a reference frame 
moving on the orbital path. In reality, attitude and orbital dynamics are coupled 
and environmental torques produced by gravity gradient, aerodynamic, and solar 
radiation pressure depend on the spacecraft’s orientation. [Ref. 6: p. 7] 



A. DESIGN PHILOSOPHY 

The goal of control system design is to cause the output variable of some 
dynamic system or process to follow a desired reference variable accurately in 
spite of changes in this reference variable, the external disturbances applied, and 
any changes in the dynamics of the process itself. Prior to beginning any control 
system design, a mathematical model of the system to be controlled is 
constructed. [Ref. 7: p. 17] The dynamic process to be controlled in this thesis 
is the attitude of the CER plus any acquired target and its equations of motion 
have been developed in the previous chapter. 

The process of regulation defines a situation in which the output variable of 
some dynamic process must follow a constant, usually zero, reference variable. 
[Ref. 7: p. 107]. The attitude control of the CER can be defined as a regulation 
process. While maneuvering or after acquiring a target, the CER's attitude must 
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remain constant despite the presence of external disturbance forces. After target 
acquisition, the CER must perform a reorientation before initiating a 
translational maneuver that will return it to the Space Station. In this 
reorientation maneuver, the reference variable is the desired orientation and the 
output variable is the current orientation. A regulation process is needed to 
reduce the orientation error, which represents the difference between the current 
and desired orientation, to zero in a timely fashion. 

This thesis designs closed-loop or feedback control systems because of then- 
inherent ability to reject disturbances and errors in the model of the dynamic 
process to be controlled. This decision, by definition, requires the introduction 
of an output sensor which can introduce noise into the control system [Ref. 7: 
pp. 107-113]. The availability of quality sensors to measure angular position and 
angular velocity is assumed in this research. Wertz [Ref. 4: pp. 155-201] 
describes in detail the various types of hardware available to determine a 
spacecraft's angular position and angular velocity. The most common 
instruments used are the rate gyro and rate-integrating gyro. Since both 
instruments have a long history of operation, and are relatively inexpensive, the 
assumption that quality sensors exist is, therefore, very reasonable. 

B. SPECIFICATIONS AND CONTROL DEVICES 

The typical feedback control system is designed such that the overall dynamic 
system response meets a set of predetermined specifications. These 
specifications, more often than not, must be translated into terms more readily 
understood by the control engineer. 

One set of specifications is known as time domain specifications and include 
such information as settling time, maximum overshoot, and damping ratio. In 
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terms of the CER’s dynamics, this refers to how long it takes for a reorientation 
maneuver and how much orientation angle overshoot is allowed before settling 
down to the final desired value. Typical spacecraft attitude reorientation control 
requires that no overshoot occurs and this corresponds to a damping ratio of 1 .0. 
Reference 1 does not specifically address the issue of settling time for CER 
attitude reorientation maneuvers or attitude regulation in the presence of 
disturbances. However, MDAC [Ref. 2: p. 9] defines a total deployment of time 
of 120 minutes. This is certainly an upper limit for any reorientation 
maneuvers. A more realistic figure is obtained from MDAC's [Ref. 2: p. 17] 
definition of major mission phases as a function of time. In this mission phase 
sequence, ten minutes is allocated to target capture. In this thesis, 70 seconds was 
selected as a reasonable settling time to accomplish any reorientation maneuver. 
Any attitude regulation, after the CER is subject to a disturbance, must be 
accomplished in a fraction of this time. A more detailed discussion of 
specifications and performance will be given in the following sections on attitude 
hold and slewing maneuvers. 

In order to achieve control over the CER some type of control device must 
be selected. Wertz [Ref. 4: p. 201-210] discusses in detail the types of devices 
available for spacecraft control. The most widely used type of device are gas jets 
or thrusters and these are the control devices that MDAC selected in their 
baseline design of the CER. 

The choice of the control actuators is, furthermore, tied to the desired 
specifications and to what is practically available. MDAC, for example, has 
chosen 1 lbf cold nitrogen thrusters. Wertz [Ref. 4: p. 206] mentions that gas 
jets are classified as cold gas and hot gas. Hot gas jets typically produce high 



25 



thrust levels (>5 N or 1.12 lbf) but rely upon a chemical reaction which must 
reach steady state. Cold gas systems produce lower thrust levels (< 1 N or 0.225 
lbf) but do not rely upon a chemical reaction which must reach steady state. 
Cold gas systems, therefore, provide more precise control and can be used 
effectively in a pulsed mode. This thesis assumes the use of the cold gas thrusters 
selected by MDAC. These thrusters are commercially available and can operate 
in a pulsed mode. 

C. ATTITUDE HOLD 

1. Optimal Control Theory 

A linear feedback control law is defined in the following form: 



u = — KX 



(3.1) 



where: 

1 . u is the applied control effort; 

2. K is a gain matrix that must be determined; 

3. X is the vector containing all the state variables, assuming that they are all 

available by either measurement or estimation. 

The gain matrix K can be determined by choosing appropriate Laplace domain 
closed-loop pole locations based on some given time domain specifications. This 
method, known as pole placement , only works well for single-input-single-output 
(SISO) systems. In multi-input-multi-output (MIMO) systems, such as the CER, 
this technique does not lead to the development of a unique control law. In 
optimal control theory , the gain matrix K is determined by minimizing a 

specified performance criterion or cost function. [Ref. 8: pp. 337-338] 
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2. Linear Quadratic Regulator 

The linear quadratic regulator is an optimal control law of the form 
shown in equation (3.1). It is called linear because the control law is a linear 
function of the system states and it is called a regulator because this type of 
control law is well suited for regulation type problems. The quadratic 
description relates to the fact that the gain matrix K, of equation (3.1), is 
calculated by minimizing a quadratic integral cost function. For the continuous 
state space system described by equation (2.23), a quadratic integral cost function 
can be formulated as: 



J = [X T (x)QX(x) + u T (x)Ru(T)]dx (3.2) 

where Q and R are symmetric weighting matrices that must be chosen by the 
control system designer. Q penalizes deviation of the state vector X from the 
origin and R penalizes the use of too much control effort. [Ref. 8: pp. 339-340] 
During attitude hold of the CER, the goal is to reject all disturbances 
and maintain a constant attitude with zero angular velocity. If all the states are 
considered initialized at zero, then the control effort must drive all states towards 

the origin after a disturbance causes a deviation from this situation. A typical 
value for Q that will cause the position angles to go to the origin and yet limit 

the angular velocities is 
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(3.3) 



Q= 



'1 0 0 0 0 0 " 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 c 0 0 

0 0 0 0 c 0 

0 0 0 0 0 c 



where c is a number less than one [Ref. 8: p. 340]. 

The choice of a proper control weighting matrix R also requires 
careful consideration. Unless a cost or penalty is imposed for using too much 
control, the design that emerges may generate control signals that can not be 
achieved by the actuators or control devices. The resulting control signal then 
saturates at the maximum signal value that can be produced and this produces, in 
most cases, the fastest possible response. The fastest possible response may be 
highly desirable but the closed-loop behavior of a system in saturation may be 
quite different from the closed loop behavior predicted by a system not in 
saturation. The system may even become unstable when the control system 
saturates and because of this consequence, R should be chosen to avoid 

saturation. [Ref. 8: p. 341] 

Hansen [Ref. 3] designed a control law based on saturating the CER's 
control input. This control scheme produced a minimum time solution but it 
required the control thrusters to be either on or off; this is known in the 
literature as bang-bang control. Some of the worst case target scenarios 
produced unstable results and this caused Hansen to choose larger values and 
different locations for the CER thrusters in order to avoid an unstable situation. 

The gain matrix K, that minimizes equation (3.2), is found by solving 

the Riccati Equation. It is, in general, a time varying matrix that, given fairly 
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general conditions which apply to the CER, eventually reaches a steady state 
value. This optimum gain matrix in the steady state is given by: 

K SS = R _1 B T M (3.4) 

where M satisfies the algebraic equation also known as the algebraic Riccati 
equation : 



0 = M A + A t M- M BR' 1 B t M + Q. (3.5) 

[Ref. 8: pp. 345-346] Many software packages, including the program 

MATLAB used in this thesis, contain subroutines that calculate the steady state 
value of K for a given dynamic system and cost function. 

Attitude hold for the CER is, therefore, accomplished by using a Linear 
Quadratic Regulator to drive the states of the system to the origin after 
experiencing some external disturbances. The gain matrix K is determined by 

supplying a MATLAB subroutine with a model of the CER’s dynamics and 
appropriately chosen Q and R weighting matrices. The steady state value of K 

is used for simplicity since K only varies near the final time and this situation is 
not encountered during most of the CER's mission. 

D. SLEWING MANEUVERS 

The problem of reorienting a spacecraft from one rest orientation to another 
rest orientation, although a problem of regulation , requires the formulation of a 
quite different control law. The linearized equations of motion, used for attitude 
hold design, are no longer valid. Slewing over a potentially large range of 
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orientation angles requires the use of the complete nonlinear equations of motion 
defined by equations (2.13) and (2.16). A linear control law may not be 
adequate for this task and the formulation of a nonlinear control law may be 
required. In addition, any slewing maneuver should, ideally, be an optimal 
maneuver. Optimal in this context refers to taking the shortest angular path. 

1. Eigenaxis Rotations 

Many spacecraft attitude control systems are currently based on a 
sequence of rotational maneuvers about each control axis. This is a natural bias 
based on the popularity of Euler angles for describing rigid body orientation. 
However, the maneuver time of such successive rotations is two to three times 
longer than that of a single maneuver about the eigenaxis. This eigenaxis is the 
principal axis developed in Euler's Principal Rotation Theorem. Euler, 
moreover, proved that the principal angle <j) of equation (2.11) is always smaller 
than the algebraic sum of three successive Euler angles and represents the 
shortest angular path between two relative orientations. Therefore, a control law 
that causes a spacecraft to reorient itself by rotating about the eigenaxis will be 
executing an optimal maneuver. [Ref. 9: pp. 375-376] 

2. Quaternion Feedback Regulator 

Reference 9 develops a nonlinear control law that takes into account the 
complete nonlinear attitude equations of motion and executes an eigenaxis 
rotational maneuver. This development begins by defining a quaternion error 
which represents the attitude error between the current orientation and the 
desired or commanded orientation: 
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Poe Poc Pic ? 2 c P3c Po 

Pic _ “Pic Poc P3c ~ P 2 C Pi 

P2c ”$2c “ P3c Poc Pic P 2 

_P3c. P 3 C ?2c “Pic Poc _ _P3_ 



(3.6) 



where terms in the matrix represent the commanded orientation expressed in 
Euler parameters. A three-dimensional error vector can be formed by 
extracting the three vector components of equation (3.6) as: 



where the first term is a nonlinear body-rate feedback term that counteracts the 
gyroscopic coupling torque found in equation (2.16), the second term is a linear 
body-rate feedback term, the third term is a linear error-quaternion or Euler 
parameter-error feedback term, and D and K are 3x3 constant gain matrices to 

be properly determined. [Ref. 9: p. 376] 

To complete the control law in equation (3.8), the matrices D and K 

must be determined. Reference 9 considered the gain selections 




(3.7) 



The complete quaternion regulator feedback control law is: 



u = colco - Deo - KP c 



(3.8) 



K = kl; 



(3.9) 
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D = dl 



(3.10) 



where k and d are positive scalars, and I is the spacecraft moment of inertia 
tensor. Global stability via Lyapunov stability analysis was proved for the 
control law of equation (3.8) provided that 

KlD>0. (3.11) 

Equation (3.11) is always guaranteed with the gains defined by equations (3.9- 
3.10). 

With global stability guaranteed for the control law in equation (3.8), 
all that remains is a proper choice of k and d. Let X be a unit vector along the 
eigenaxis. The orientation is then expressed as: 

g = sin(<t>/2)L (3.12) 

Substituting equation (3.8) into equation (2.16) yields the following closed-loop 
equation: 



cb = -I 1 wlco + 1 1 (wlco - Deo - K p) . 



(3.13) 



Further substitution of equations (3.9-3.10) produces: 

co = r‘(-dIco-kIp) (3.14) 
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or by rearranging terms: 



Icb + dlco + kip = 0 (3.15) 

When the angular rate co is small, and when an eigenaxis rotation is assumed, the 
angular rate can be approximated as: 



co = 4 )i. (3.16) 

Further substitution of equations (3.16) and (3.12) into equation (3.15) yields: 

($+d<j> + ksin(<t>/2))lX. = 0. (3.17) 

Since the moment of inertia I is, by definition, a positive definite matrix and the 
unit vector X is non zero, then ]X * 0 and equation (3.17) becomes: 

6 + d<t> + ksin((j>/2) = 0. (3.18) 

If the sin(<j>/2) is approximated by (j)/2, equation (3.18) is further reduced to: 

6 + d<j) + k<j)/2 = 0. (3.19) 

Equation (3.19) is the well-known linear second order equation where the 
damping ratio £ and the natural frequency co N satisfy: 
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d = 2Cco N ; 



(3.20) 



k = 2co*. (3.21) 

Therefore, proper selection of the damping ratio and the natural frequency 
defines the positive scalars d and k. [Ref. 9: p. 377-378] 

Since this thesis has assumed a damping ratio of one, a required 
maneuver settling time is converted to a required natural frequency by: 

T s =4/Ccd n . (3.22) 

To account for the nonlinear effect of sin(<j)/2) when <j) is large, equation (3.22) 
was modified by Ref. 9 as: 



T s = 8/Cco n . (3.23) 

Therefore, a quaternion feedback regulator control law is defined by a desired 
orientation and by a desired maneuver settling time as: 



d = 2 co n = 16/ T s ; 


(3.24) 


k = 2a>N =128/T S 2 ; 


(3.25) 


u = colco- (16 / T s )Ico - (128 / T s 2 )IP„ . 


(3.26) 
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IV. ADAPTATION LAW DESIGN 



Knowledge of the moment of inertia tensor for the CER and acquired target 
is required by both the control laws developed in Chapter III. In the case of the 
large angle control law, imprecise knowledge of the moment of inertia tensor 
results in a non-eigenaxis slewing maneuver. An argument can also be made that 
precise knowledge of the moment of inertia tensor should limit the amount of 
thruster firings required and save propellant fuel. Therefore, the ability to 
estimate the moment of inertia tensor and provide this information to the control 
algorithms is a very beneficial addition to the previously developed control 
algorithms. 

A. ADAPTIVE CONTROL THEORY 

An adaptive controller differs from a static or ordinary controller in only 
one respect. In an adaptive controller, the controller parameters are variable 
and there is a mechanism for adjusting these parameters on-line based on signals 
available from the overall system. The two main approaches for constructing 
adaptive controllers are: the model reference method and the self -tuning 

method. A schematic representation of each of these methods is presented in 
Figure 7. [Ref. 10: p. 315] 

This thesis employs the self-tuning method. The design of an adaptive 
controller by the self-tuning method involves choosing a control law based on 
variable parameters and choosing an adaptation law for adjusting those 
parameters. The controller, therefore, couples a previously designed control law 
with an on-line parameter estimator. [Ref. 10: pp. 319-323] The previous 
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Figure 7. Adaptive Control Methods [Ref. 10: pp. 315, 320] 
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chapter developed two control laws based on knowledge of the CER moment of 
inertia tensor; the next section develops an adaptation law based on the Kalman 
filter equations. 

B. KALMAN FILTER DESIGN 

1. General Kalman Filter Equations 

The Kalman filter is an observer or state estimator for a dynamic 
process given by the following discrete state space equation: 

X(k + 1) = OX(k) + Aj u(k) + A 2 w(k) (4.1) 



where: 

1. X(k) is the state vector at the present time k; 

2. X(k + 1) is the state vector one time step in the future; 

3. O is the discrete-time version of the A matrix given in equation (2.23); 

4. u(k) is the applied control; 

5. w(k) is an unknown random input called the plant driving noise\ 

6. Aj is the discrete version of the B matrix given in equation (2.23); 

7. A 2 is the random input influence matrix and is often identical to A, . 

In addition, the measurements of the system are given by: 

Y(k) = CX(k) + v(k) (4.2) 

where v(k) is a random vector known as measurement noise. [Ref. 11: p. 27] 
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The Kalman filter provides an optimal estimate of the system described 
by equation (4.1) and (4.2) by minimizing the mean square error between the 
actual states and the estimated states. The Kalman filter has the following form: 

X(k + Ilk + 1) = X(k + Ilk) + G(k + l)[Y(k + 1) - Y(k + ilk)] (4.3) 



where: 

1 . the notation (k + ilk) refers to the discrete value at time k + 1 based on data 
accumulated through time k; 

2. X(k + Ilk) is the estimate of the states given data through time k; 

3. Y(k + Ilk) is the estimate of the measurements given data through time k; 

4. G(k + 1) is known as the Kalman filter gain; 

5. X(k + Ilk + 1) is estimate of the states given data through time k + 1. 

A set of recursive equations, in the proper order, that solve for the Kalman filter 
gain and equation (4.3) are: 



P(k + Ilk) = OP(klk)O x + A 2 WA 2 t ; 


(4.4) 


G(k + 1) = P(k + llk)C T [CE(k + ll k)C T + v]~' ; 


(4.5) 


P(k + ll k + 1) = [l - G(k + l)C]P(k + ll k) ; 


(4.6) 


X(k + ll k) = OX(kl k) + A, u(k) ; 


(4.7) 


Y(k + ll k) = CX(k + ll k) ; 


(4.8) 
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X(k + Ilk + 1) = X(k + Ilk) + G(k + l)[Y(k + 1) - Y(k + Ilk)]; 



(4.9) 



where: 

1 . P is the covariance matrix of the estimation error; 

2. W is the covariance matrix of the plant driving noise; 

3. V is the covariance matrix of the measurement noise; 

4. 1 is an appropriately dimensioned identity matrix. 

Note that final implementation of equations (4.4-4.9) requires an initial estimate 
of the states X(OlO) and an initial estimate of the covariance matrix P(OlO) [Ref. 

11: pp. 27-29]. The initialization of equations (4.4-4.9) is discussed in Chapter 
V. 

2. Linear Model 

For attitude hold, the linear quadratic regulator requires a model of the 
CER's dynamics and this is provided by equations (2.23-2.25). The B matrix 
defined by equation (2.24) contains the elements of the CER's inverse moment of 
inertia tensor which is the system parameter that needs to be estimated. 

To apply the previously defined Kalman filter equations, a plant 
equation similar to equation (4.1) is formulated as: 

I" 1 (k + 1) = I" 1 (k) + w(k) (4.10) 

where I _1 (k) is a vector that contains the six independent elements of the inverse 
moment of inertia tensor and w(k) is random plant driving noise. By comparing 
equations (4.1) and (4.10), the following correspondences are noted: 
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(4.11) 



0 = 1 ; 

A, =0; 
A, = l. 



(4.12) 

(4.13) 



This basically assumes that the inverse moment of inertia does not change, other 
than by the addition of some random noise, from one time step until the next. 

An accompanying measurement equation is formulated from a 
linearized version of equation (2.16) as 

co = r‘u. (4.14) 



Equation (4.14) can be algebraically rearranged into the following form: 



Ui (k) 



w(k) = 



0 

0 



u 2 (k) 

u,(k) 

0 



u 3 (k) 0 

0 u 2 (k) 

Ul (k) 0 



0 

u 3 (k) 

u 2 (k) 



0 

0 



I- 1 (k) + y(k) 



u 3 (k) 



(4.15) 



where v(k) is the combined measurement noise that results when angular 
acceleration and applied control torques are measured. By comparing equations 
(4.2) and (4.15), the following correspondence is noted: 



C(k) = 



u,(k) 

0 

0 



u 2 (k) 

u,(k) 

0 



u 3 (k) 0 

0 u 2 (k) 

u 3 (k) 0 



0 

u 3 (k) 

u 2 (k) 



0 

0 . 
u 3 (k) 



(4.16) 
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The previously defined Kalman filter can now be used to estimate the 
inverse moment of inertia tensor of the CER plus any acquired target. This on- 
line estimation assumes that: 

1. Reasonable angular acceleration measurements or estimates are available; 

2. The applied control torques are available; 

3. The covariance matrices for plant and measurement noise can be computed; 

4. Initial estimates for the inverse moment of inertia and covariance matrix of 
estimation error can be computed. 

All the above assumptions are discussed further in Chapter V. 

3. Nonlinear Model 

An estimate of the moment of inertia for the CER and any acquired 
target is required by the quaternion feedback regulator in order to accomplish 
any slewing maneuvers. A plant equation similar to equation (4.10) is formed 
as: 



where equations (4.11-4.13) are still valid. Since large angle motion is 
inherently nonlinear, it is only prudent to form a measurement equation from the 
original nonlinear version of equation (2.16) which has been rearranged as: 



I(k + 1) = I(k) + w(k) 



(4.17) 



U = C0lG) + Icb. 



(4.18) 



Further algebraic manipulation results in the following form: 



u(k) = C(k)I(k) 



(4.19) 
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where: 



cbi 


C0 2 -0)^3 


cb 3 + C0 1 C0 2 


-0)30)2 


2 2 
0) 2 -0)3 


0) 3 0) 2 




CO1CO3 


COi +C 0 2 Gl >3 


2 2 
CO3 -cof 


Cb 2 


0)3 -0)j0) 2 


-0)i0) 3 


( 4 . 20 ) 


-co 1 co 2 


cof -co 2 


COj -0) 2 0) 3 


0)j0) 2 


d) 2 + 0) 1 0) 3 


0) 3 





and the time dependence of each term is assumed. Note that while the dynamics 
of the CER are nonlinear, the dynamics of the CER's moment of inertia remain 
linear. 

The previously defined Kalman filter can now be used with this model 
to estimate the moment of inertia tensor of the CER plus any acquired target. 
This on-line estimation scheme utilizes the same assumptions listed in the 
previous section. 
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V. RESULTS 



A. SIMULATION PROGRAM AND IMPLEMENTATION 

Attitude control performance for the CER and target was evaluated through 
computer simulation. The simulation programs were written in MATLAB and 
are listed in the Appendix. The continuous-time equations of motion given in 
Chapter II and the continuous-time control algorithms given in Chapter III were 
simulated using discrete-time versions of these equations. A small sampling 
period of 75 milliseconds was chosen to ensure that the discrete equations 
accurately represented the continuous-time equations. 

The selection of thruster size and location inherently limits the amount of 
available control torques. To ensure realistic simulation results, hard limits have 
been included in the simulation programs so that control signals greater than 
those given in Chapter I are not generated. 

Plant noise has been added into the simulation programs to further increase 
the realism of the computer simulation results. Typical spacecraft disturbance 
torques are due to thruster misalignment and solar pressure. Kaplan [Ref. 12: p. 
241] lists some assumed values for disturbance torques. The maximum 
magnitudes of these assumed values are: 

1. Thruster misalignment torque: 8.5 x 10~5 N-m or 6.27 x 10~5 ft-lbf; 

2. Solar Pressure torque: + 1.0 x 10"4 N-m or + 7.4 x 10"5 ft-lbf. 

As a worst case guess, plant noise was modeled as a random three dimensional 
vector with a normal or Gaussian distribution and a standard deviation equal to 
the sum of these two disturbance torques, 1.4 x 10'4 ft-lbf. 
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The control laws developed in Chapter II assume the availability of sensors 
to measure angular position and angular velocity. The sensors have limited 
accuracy and introduce measurement noise into the control system. Wertz [Ref. 
4: pp. 199-200] states that typical rate gyros have a resolution of less than 0.01 
degree per second and typical rate-integrating gyros have a resolution of 0.003 
degree. Reference 13 indicated that accuracies of 0.01 degree and 0.0003 degree 
per second are quite reasonable. As a worst case guess, measurement noise was 
modeled as two random three dimensional vectors with a normal or Gaussian 
distribution and standard deviations equal to 0.003 degree and 0.0001 degree per 
second. 

The linear quadratic regulator used for attitude hold required a choice of Q 
and R weighting matrices. After some trial and error, weighting matrices were 
obtained that allowed the states to be driven towards the origin in a reasonable 
amount of time (ten seconds) and also avoided saturation of the applied control 
torques: 



T 

0 




0 

0 



0 0 0 
1 0 0 
0 1 0 
0 0 0 
0 0 0 
0 0 0 



0 0 “ 
0 0 
0 0 
0 0 ’ 
0 0 
0 0 



R = 



2x10“ 

0 

0 



2x10“ 



0 

0 

2xl0" 5 



(5.1) 



(5.2) 
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Finally, the adaptation algorithm developed in Chapter IV required some 
computations and assumptions. Reference 13 indicated that angular 
accelerometers exist for spacecraft applications but are very expensive and not as 
common as rate gyros and rate-integrating gyros. In this thesis, angular 
acceleration measurements were generated by numerically differentiating the 
given angular velocity measurements as follows: 

o)(k) = (co(k + l)-o)(k))/T i (5.3) 

where T g is the sample period or interval. Alternative methods for generating 
angular acceleration measurements might involve filtering the already available 
angular velocity measurements or using an extended Kalman filter to estimate the 
moment of inertia tensor and the angular accelerations simultaneously. The 
initial covariance of estimation error matrix P(OlO) was calculated empirically by 

computing the total moment of inertia tensor for several targets in various 
locations within the CER's net and then computing the overall mean and standard 
deviation. The covariance matrix of plant driving noise W and the covariance 
matrix of measurement noise V were computed as follows: 



W = wl; 


(5.4) 


V = vl; 


(5.5) 



where 1 is an appropriately dimensioned identity matrix, and v and w are scalar 
values equal to the variances of the assumed plant noise and measurement noise. 
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The plant is the CER and target moment of inertia and this author assumes that it 
does not change very significantly. A small number should, therefore, be chosen 
for w, with zero an ideal choice. However, simulations indicated that a zero 
plant noise covariance matrix leads to unsatisfactory results. In this thesis, the 
following small number was used: 



w = 1x10 8 . (5.6) 

The measurement noise takes into account the errors and noise that are 
introduced by thruster misalignment, and measuring angular acceleration and 
angular velocity. After some trial and error via computer simulations, the value 
chosen empirically for v was: 



v = lxlCT 4 . (5.7) 

The adaptation scheme constantly updates the control algorithms as it 
estimates the moment of inertia. The Kalman filter has some inevitable 
transients prior to convergence and computer simulation results indicated that the 
Kalman filter may even provide an estimate for the moment of inertia that is not 
physically realistic. Since the moment of inertia tensor is by definition a positive 
definite matrix (the eigenvalues are always positive for a real physical body), a 
test has been incorporated within the control and adaptation schemes: estimates 
of the moment of inertia are only passed to the control law if all the eigenvalues 
are positive. 
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Hansen [Ref. 3: p. 21] computed the total moment of inertia tensor for the 
CER and several worst case target capture scenarios. These moments of inertia, 
tabulated in Table 1, are used in the following sections to evaluate the control 
laws developed in this thesis. 



TABLE 1. TARGET CAPTURE MOMENT OF INERTIA TENSORS 



Case/Description 


Moment of Inertia Tensor 
(slug-foot 2 ) 


1. No Target 




"39.6 0 0" 

0 55 0 

0 0 55 




2. Man + manned maneuvering 
unit in net center not aligned with 
CER. 




"112.9 2.4 -111.9" 

2.4 534.9 6.4 

-111.9 6.4 497.6 




3. 850 lb point mass at net edge 
along the x-axis. 




' 69.3 0 -178.2' 

0 1153.8 0 

-178.2 0 1124.1 




4. 850 lb point mass at net edge 
along the y-axis. 




' 98.9 -110.5 -110.5' 

-110.5 496.1 -29.7 

-110.5 -29.7 496.1 




5. 850 lb point mass at net edge 
along the z-axis. 




' 369.5 0 -368.4" 

0 796.4 0 

-368.4 0 466.4 




6. 850 lb point mass at X=Z, Y=1 
from net center. 




" 172.7 -93.7 -282.3" 

-93.7 839.7 -39.8 

-282.3 -39.8 732.9 
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B. ATTITUDE HOLD 



Evaluation of attitude hold for the CER and acquired target was 
accomplished by initializing the position angles and angular velocities at some 
reasonable values and then allowing the control system to drive these values 
towards the origin. This thesis employed the same the initial conditions as 
Hansen [Ref. 3: p. 43]: 2.0 degrees for all position angles and 0.2 degrees per 
second for all angular velocities. Throughout the attitude hold simulations, the 
adaptive control system has been initialized with the moment of inertia defined 
by Case 1 in Table 1. 

1. Comparison of Adaptive and Nonadaptive Control 

The control system developed in Chapter III was compared with a 
nonadaptive version of the same control system for a CER and target moment of 
inertia defined by Case 2 in Table 1. The nonadaptive controller was given a 
constant value for the CER and target moment of inertia: defined by Case 1 in 
Table 1. The simulation results are listed below in Table 2. 



TABLE 2. ADAPTIVE/NONADAPTIVE SIMULATION RESULTS 



Controller Design 


Settling Time (sec) 


Fuel Used 


1. Nonadaptive, 
1= Case 1 


30 


1009 


2. Adaptive 


10 


398 



A figure of merit for any spacecraft control system design is the amount of fuel 
used. In this thesis, the sum of the absolute value of all control used for each 
design was calculated. This value is listed as fuel used in Table 2, although it is 
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actually only proportional to the amount of fuel used. Note that the settling time 
in Table 2 is the maximum settling time for all three position angles. 

Simulation plots for the nonadaptive design (1= Case 1) and the adaptive 
design are displayed in Figure 8. Both designs provide for an overall stable 
closed loop system. However, the nonadaptive design is clearly much more 
oscillatory, requires a long period for all oscillations to completely dampen out, 
and requires more control effort. The adaptive design starts out with the worst 
possible guess (the CER alone) for the moment of inertia but rapidly and 
correctly estimates the actual moment of inertia of the CER and target. The final 
estimate is very close to the actual moment of inertia defined by Case 2 in Table 
1 : 



‘ 112.7 


2.8 


-111.5' 


2.8 


540 


3.3 


-111.5 


3.3 


498.1 
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Figure 8. Simulation Results for Attitude Hold 
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2. Adaptive Control with various target scenarios 

The adaptive control system developed for attitude hold in Chapter III 
was simulated for all the target-capture moment of inertia tensors defined by 
Hansen [Ref. 3]: Cases 1-6 in Table 1. The simulation results for this adaptive 
attitude control design and the nonadaptive minimum time controller designed by 
Hansen [Ref. 3] are listed below in Table 3. 



TABLE 3. ATTITUDE HOLD SIMULATION RESULTS 



Target Case 


Settling Time (sec) 
Adaptive Design 


Settling Time (sec) 
Hansen Design 


1 


3.0 


1.5 


2 


10.0 


6.6 


3 


15.0 


oo 


4 


8.0 


oo 


5 


15.0 


15.6 


6 


13.0 


CO 



The adaptive control system was stable and exhibited reasonable settling times 
for all target cases. The minimum time control system designed by Hansen [Ref. 
3] exhibited unstable results for three of the target capture cases: as depicted by 
the infinite settling times. 

Table 4 lists the final estimate of the CER and target moment of inertia 
tensor for each of the target cases. The results are close but not exact. The 
differences may be attributed to the use of noisy measurement signals and the 
short time (10 seconds) available for system measurements. As the control 



51 



torque and angular acceleration measurements become very small, the estimation 
scheme is essentially provided with no new information and therefore can not 
accurately estimate the moment of inertia tensor. 



TABLE 4. MOMENT OF INERTIA ESTIMATES 



Target 

Case 


Actual Moment of 
Inertia Tensor 


Estimate of Moment 
of Inertia Tensor 


1 




"39.6 0 0" 

0 55 0 

0 0 55 






"39.6 -0.3 0.3" 

-0.3 55.2 0.2 

0.3 0.2 54. 4_ 




2 




"112.9 2.4 -111.9" 

2.4 534.9 6.4 

—111.9 6.4 497.6 






"112.7 2.8 -111.5] 

2.8 540 3.3 

-111.5 3.3 498.1 J 


3 




" 69.3 0 -178.2' 

0 1153.8 0 

-178.2 0 1124.1 






' 69.3 -7.0 -179.4" 

-7.0 1155.8 -4.4 

-179.4 -4.4 1131.6 




4 




' 98.9 -110.5 -110.5' 

-110.5 496.1 -29.7 

-110.5 -29.7 496.1 






98.7 -110.6 -109.9' 

-110.6 284.5 182.0 

-109.9 182.0 282.7 




5 




' 369.5 0 -368.4" 

0 796.4 0 

-368.4 0 466.4 _ 






"411.4 5.0 -417.0" 

5.0 797.0 -5.9 

-417.0 -5.9 522.9 _ 




6 




" 172.7 -93.7 -282.3" 

-93.7 839.7 -39.8 

-282.3 -39.8 732.9 _ 






' 173.7 -90.2 -286.1" 

-90.2 851.5 -52.6 

-286.1 -52.6 746.8 _ 





C. SLEWING MANEUVERS 

Evaluation of slewing maneuvers for the CER and acquired target was 
accomplished by: initializing the position angles and angular velocities at zero, 
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selecting a final orientation in terms of the chosen 3-2-1 Euler Angles , and 
selecting a desired maneuver settling time. The adaptive quaternion feedback 
regulator developed in Chapters III and IV was compared with two nonadaptive 
quaternion feedback regulators that assumed a CER and target moment of inertia 
as defined by Case 1 in Table 1. The adaptive quaternion feedback regulator was 
initialized with this same moment of inertia. A desired orientation of 50 degrees 
for each position angle was selected for each controller design. The simulation 
results for a target defined as Case 2 are listed below in Table 5. 



TABLE 5. SLEWING MANEUVER SIMULATION RESULTS 



Controller 


Desired 


Desired 


Actual 


Max 


Fuel 




Settling 


Damping 


Settling 


% 


Used 




Time (sec) 


Ratio 


Time (sec) 


Overshoot 




Adaptive 


70 


1.0 


70 


0 


826 


Non- 


70 


1.0 


+100 


55 


722 


adaptive 












Non- 


11 


2.5 


70 


4 


958 


adaptive 













The adaptive design clearly meets the desired settling time and desired 
overshoot requirements. The first nonadaptive design appears to use more fuel 
than the adaptive design, but it is very oscillatory and does not meet the desired 
settling time. A second nonadaptive design was generated by adjusting the 
desired settling time and damping ratio, which changes the parameters k and d in 
the control law developed in Chapter III, until the actual settling time approached 
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the original desired settling time of 70 seconds. This second nonadaptive design 
meets the required settling time but the response is more oscillatory and uses 
more fuel: as compared with the adaptive design. Figures 9, 10, and 11 further 
clarify the differences between the adaptive and nonadaptive designs. 

The simulation plots of the quaternions or Euler parameters can be used to 
check for an eigenaxis rotation: which is defined by Ref. 9 as a straight line in 
any plot of the quaternions or Euler parameters. The first nonadaptive 
controller is clearly executing a noneigenaxis rotation. The second nonadaptive 
design controller and the adaptive controller are, almost exactly, executing an 
eigenaxis rotation. The adaptive controller is, therefore, the only design that: 

1 . minimizes the desired position angle overshoot; 

2. meets the desired settling time; 

3. minimizes the fuel used; 

4. executes an eigenaxis rotation; 

5. has the ability to react to different and changing moments of inertia. 
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Figure 9. Simulation Results for Nonadaptive (Design 1) Slewing 
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Figure 10. Simulation Results for Adaptive Slewing 
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Figure 11. Simulation Results for Nonadaptive (Design 2) Slewing 
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VI. CONCLUSIONS 



A. ATTITUDE HOLD 

The adaptive attitude control system developed in this thesis to maintain 
attitude hold for the CER and an acquired target provided superior control 
during the capture of all previously defined worst case target scenarios. This 
adaptive linear quadratic regulator provided stable results with very reasonable 
settling times, and used a modest amount of fuel as compared with the previously 
designed nonadaptive minimum time control system [Ref. 3] and a nonadaptive 
linear quadratic regulator design. The additional computational burden of 
adaptive control is, therefore, compensated by a substantial improvement in 
control system performance. 

B. SLEWING MANEUVERS 

Fairly accurate knowledge of the moment of inertia tensor is essential for 
slewing maneuvers that are well damped, accomplished in a timely manner, use 
minimal fuel, and are executed as eigenaxis rotations. The adaptive quaternion 
feedback regulator developed in this thesis clearly provides this type of slewing 
maneuver control and outperforms nonadaptive quaternion feedback regulators. 
Again, the additional computational burden of adaptive control is more than 
compensated for by a substantial improvement in control system performance. 

C. FUTURE RESEARCH 

Both adaptive control system designs developed in this thesis are practical 
designs that will result in a more reliable and fuel efficient Crew/Equipment 
Retriever for the Space Station Freedom. These designs, however, can also be 
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applied to more general spacecraft attitude control problems. The adaptation 
scheme can easily be altered to account for spacecraft moments of inertia that are 
not only unknown but are subject to frequent and significant changes. In this 
respect, the adaptive linear quadratic regulator is an ideal candidate for large 
space structures of the future. Instead of limiting space structure design to 
shapes and load distributions that lead to quick approximations of the overall 
moment of inertia and minimize any changes, this adaptive control system would 
remove all these restrictions and become, in current science fiction terminology, 
an automatic attitude damping system. 

Future research could examine alternative adaptation schemes and employ 
filtering techniques to obtain angular acceleration measurements from angular 
velocity measurements, since the latter are more generally available. Additional 
computer simulations could attempt to model the CER and target as a flexible 
structure and account for any target movement within the CER’s capture net 
mechanism. 
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APPENDIX. MATLAB SIMULATION PROGRAMS 



% PROGRAM NAME: CERSimNLl 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: actualmoi, CERLQR, guessmoi 

% 

% DESCRIPTION: 

% The below program will simulate the CER Control System for Small 
% Angle Motion or Attitude Hold. The Control Law is a Steady State Linear 
% Quadratic Regulator. However, it is an Adaptive Controller in that the 
% Inverse MOI of the CER + Target is unknown and potentially variable. A 
% KALMAN Filter is used to estimate the Inverse MOI and then update the 
% feedback gains produced by the Linear Quadratic Regulator (LQR). It is 
% assumed that the angular accelerations are available as measurements. 

% 

clear 

% Prompt user for first guess inverse MOI tensor 
guessmoi; 

% 

% Define the initial state estimate for the Kalman Filter 
est_x=[Iig(l,l);Iig(l,2);Iig(l,3);Iig(2,2);Iig(2,3);Iig(3,3)]; 

% 

% Define the Q and r parameters for the LQR 

Q=diag([ 1,1,1 *0,0,0]); % The Q matrix weights the system states 

% 

r=2.0e-5; % This is used to weight the control inputs used 
% This is the for used R=r*eye(3,3); 

% 

% Define the actual inverse MOI tensor, as Ii. 

% This is placed into the MATLAB workspace by a separate 

% subroutine. 

actualmoi; 

% 

T=75e-3; % The sampling time. 

B=[zeros(3,3);Ii]; % The B matrix. 

% 

t(l)=l; %Initial sampling index 

% Initial conditions in degrees & degrees/sec and in radians & radians/sec 
x(:,l)=[2;2;2;0.2;0.2;0.2].*(pi/180) 

% 

y(:,l)=x(:,l); % Initial measurements 

% 
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% Initial Kalman Filter Parameters 
phik=eye(6); % Define the phi matrix for the plant 
% 

Ew=(le-8)*eye(6); % Plant Noise Covariance Matrix 
V=1.0e-4; 

Ev=V*eye(3); % Measurement Noise Covariance Matrix 
% Prompt user for Covariance Estimation Error 

ques3=input('Specify the Covariance Estimation Error, P(0/0). < 0 > The 

Default. < 1 > The Identity Matrix. < 2 > A matrix of zeros.»>’); 

% 

if ques3==0, 

oldP=[-0.01 14 -1.1 4e-6 2.92e-5 2.166e-5 4.84e-7 2.95e-5 

-1.14e-6 0.0001 -2.56e-7 -1.9e-7 -4.25e-9 -2.587e-7 

2.92e-5 -2.56e-7 -0.0026 4.864e-6 1.087e-7 6.623e-6 

le-6 le-6 le-6 -.0019 le-6 le-6 

le-6 le-6 le-6 le-6 -4.246e-5 le-6 

le-6 le-6 le-6 le-6 le-6 -2.587e-3]; 

elseif ques3==l, 

oldP=eye(6); 

else, 

oldP=zeros(6,6); 

end 

% This is P(0/0), the covariance matrix of the estimation error 

% 

% Prompt user for number of simulation steps 

ques4=input('Enter the number of Simulation Steps. < 0 > The Default, which is 
160. Otherwise, enter your own value (within reason. ...unless you have all 
day! !»>'); 
if ques4==0, 

NUM=160; 

else, 

NUM=ques4; 

end 

% 

rand(’normar); 

% 

rd=(l-4e-4)*rand(3,NUM); % Random Plant Disturbance 
mnl=(1.745e-4)*rand(3,NUM); % Random Measurement Noise 
mn2=(5e-6)*rand(3,NUM);% 
mn=[mnl;mn2]; 

% 

% Now, simulate the CER + Target 
for N=1:NUM 
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% 

wskew=[0 -x(6,N) x(5,N);x(6,N) 0 -x(4,N);-x(5,N) x(4,N) 0]; 

% Skew symmetric matrix 
a=-l *Ii*wskew*inv(Ii); 
i5=tan(x(2,N))*sin(x(l,N)); 
i6=cos(x( 1 ,N))*tan(x(2,N)); 
j5=cos(x(l,N)); 
j6=-sin(x(l,N)); 
k5=sin(x(l,N))/cos(x(2,N)); 
k6=cos(x( 1 ,N))/cos(x(2,N)); 

A=[0 0 0 1 15 i6;0 0 0 0 j5 j6;0 0 0 0 k5 k6;0 0 0 a(l,l) a(l,2) a(l,3);0 0 0 a(2,l) 
a(2,2) a(2,3);0 0 0 a(3,l) a(3,2) a(3,3)]; 

% 

% Now, discretize the state space equations 
[phi,delul]=c2d(A,B,T); 

% 

% Call the CERLQR function to determine the Steady State Feedback Gain 
% Matrix 

K=CERLQR(Q,r,Iig); 

% 

t(N+l)=N+l; % Increment time index 

% 



U(:,N)=-K*y(:,N); 

% Now, apply the limits for the control torques 

if abs(U(l,N))>4 

U(l,N)=4*sign(U(l,N)); 

else 

U(1,N)=U(1,N); 

end 

if abs(U(2,N)) > 3 

U(2,N)=3*sign(U(2,N)); 

else 

U(2,N)=U(2,N); 

end 

if abs(U(3,N)) > 3 

U(3,N)=3*sign(U(3,N)); 

else 

U(3,N)=U(3,N); 

end 

% 

torque(:,N)=U(:,N)+rd(:,N); 

% The plant discrete state equations 
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x(:,N+l)=phi*x(:,N)+delul*torque(:,N); 

% Measurement Equation 

y(:,N+l )=x(:,N+l )+mn(:,N); 

% 

xd(:,N)=A*x(:,N)+B*U(:,N); 

yka(:,N)=xd(4:6,N); % This measurement equation takes the angular velocities 
% and obtains angular accelerations (like an accelerometer ). 

% 

yk(:,N)=(y(4:6,N+l)-y(4:6,N))./T; % This equation obtains angular accelerations 
% by numerical differentiation. 

% 

% Now, call the Kalman Filter in order to estimate the 
% actual Inverse MOI and update the LQR Gain Matrix. 

% 



tl=U(l,N); 

t2=U(2,N); 

t3=U(3,N); 

Ck=[tl t2 t3 0 0 0;0 tl 0 t2 t3 0;0 0 tl 0 t2 t3]; 

% 

% Below are the Kalman Filter and gain equations 
% 

% First, the Gain Equations 

newP=phik*oldP*phik' + Ew; 

G=newP*Ck'*inv(Ck*newP*Ck' + Ev); 
oldP=(eye(6) - G*Ck)*newP; 

% 

% 

est_xx=phik*est_x(: ,N); 
est_y= Ck*est_xx; 

est_x(:,N+l)= est_xx + G*(yk(:,N)-est_y); 

% 

II l=est_x(l,N); 

I12=est_x(2,N); 

I13=est_x(3,N); 

I22=est_x(4,N); 

I23=est_x(5,N); 

I33=est_x(6,N); 

IEST=[I1 1 112 II 3;I12 122 I23;I13 123 133]; 

% 

% Check the inverse MOI guess for physical reality. Eigenvalues must be > 0. 
eigcheck(l :3,N)=[eig(IEST)]; 

% 
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% Test for physical reality ( All eigenvalues of IEST >0) 

checke(N)=l *sign(eigcheck(l ,N))+1 *sign(eigcheck(2,N))+l *sign(eigcheck(3,N)) 

% 

if checke(N)==3, 

% Update the initial guess inverse MOI 

Iig=IEST; 

% Otherwise, do not update the control law 

else, 

end 

% 

end % End of Simulation Loop 

% 

% Convert the time vector to actual time 
t=t.*T; 

% create a special time vector for the control input plot 

tU=0:N-l; 

tU=tU.*T; 

% Calculate the total fuel used (actually this is proportional to it) 

FUELUSED=sum(sum(abs(U))) 

% 

% End of Program 



% PROGRAM NAME: guessmoi 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% The below program will prompt user for first guess inverse MOI tensor. 
input('Welcome to the Non-Linear Simulation of the CER and Target. Please 
choose a first guess inverse MOI tensor for the LQR. <0 > The Default. 
MAN+MMU in Net Center with CER Frame of References 1 > The CER alone. 
< 2 > Case Two from Hansen Thesis. *Any other number will use an average 
of all the inverse MOI tensors>»'); 
s=ans; 
if s==2, 

Iig= [0.0 1 14 -0.0001 0.00256;-0.0001 0.0019 4.246e-5;0.00256 4.246e-5 
0.002587]; 

% 

elseif s==l, 

Iig=inv([39.6 0 0;0 55 0;0 0 55]); 
elseif S— — 0 

Iig=[0.0138 0 0.003;0 0.0018 0;0.003 0 0.0026] 
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% 

else, 

Iig=inv([l 43.8 -33.6 -175.2;-33.6 646 -10.5;-175.2 -10.5 562]) 
end 

% 

Ig=inv(Iig) 

% This is the guess MOI. 

% 

% End of Program 



% PROGRAM NAME: actualmoi 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% This program will question the user for a CER + target 
% moment of inertia or its inverse. 

% 

caseN=input('Define the actual Moment of Inertia by case number»'); 
if caseN ~ 1 

MOI=[39.6 0 0;0 55 0;0 0 55]; 

% This is the CER alone. 

% 

PKPlf — / 

MOI=[112.9 2.4 -1 1 1.9;2.4 534.9 6.4;-111.9 6.4 497.6]; 

% MAN + MMU 

% 

MOI=[69.3 0 -178.2;0 1153.8 0;-178.2 0 1 124.1]; 

% 

elseif caseN==4, 

MOI=[98.9 -110.5 -1 10.5;-1 10.5 496.1 -29.7;-110.5 -29.7 496.1]; 

% 

elseif caseN==5, 

MOI=[369.5 0 -368.4;0 796.4 0;-368.4 0 466.4]; 

% 

elseif caseN==6, 

MOI=[172.7 -93.7 -282.3;-93.7 839.7 -39.8;-282.3 -39.8 732.9]; 

% 

else, 

end 

Ii=inv(MOI); 

% 
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End of Program 



% PROGRAM NAME: CERLQR 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: dlqr 

% 

% DESCRIPTION: 

% The below subroutine will compute the Steady State Linear 
% Quadratic Regulator Control Gains for the CER using a Linear 
% State Space Model. It will receive updates on the inverse 
% Moment of Inertia (MOI) Tensor from another subroutine 
% that will estimate this matrix using a Kalman Filter. 

% 

% function K= CERLQR(Q,r,Ii) 

% Here are the input parameters: 

% Q is the weighting matrix for the states 
% r is the scalar used to weight the control inputs 
% Ii is the inverse MOI Tensor 
% 

function K=CERLQR(Q,r,Ii); 

% 

R=r*eye(3,3); 

% 

% Now, define the state space equations and discretize. 

T=75e~3; % The sampling time. 

A=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1;0 0 0 0 0 0;0 0 0 0 0 0;0 0 0 0 0 0]; 
% The above is the A matrix. 

% 

B=[zeros(3,3);Ii]; % The B matrix. 

% Now, discretize the state space equations 
[phi,delul]=c2d(A,B,T); 

% 

% Now, call the function dlqr to calculate the feedback gains 

% 

[K,S]=dlqr(phi,delu 1 ,Q,R); 

% END of the function 



% PROGRAM NAME: replotSA 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: subplot, plot 

% 

% DESCRIPTION: 

% Plotting Program for Small Angle Motion 
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clg; 

subplot(221); 

plot(t,x(l,:),t,x(2,:),'~',t,x(3, :),'-. '),grid 
xlabel(Time - seconds'), ylabel(’Position - radians') 
title('Position Angle Histories’) 

% 

subplot(222); 

plot(t,x(4,:),t,x(5,:),'~',t,x(6, :),'-.’), grid 

xlabel('Time - seconds’), ylabel('Angular velocity - rad/sec’) 

titlefAngular Velocity Histories') 

% 

subplot(223); 

P lot(tU,U(l,:),tU,U(2,:), , — , ,tU,U(3, :),’-.'), grid 
xlabel('Time - seconds'), ylabel(’Control Torque - foot-lbs’) 
title(’Control Torque Histories') 

% 

subplot(224); 

% Plot fuel usage plot 
plot(tU,sum(abs(U))),grid,title('Fuel Usage') 
xlabel(’Time - seconds'), ylabel('Magnitude’); 
pause; 

% Prompt user to see if additional plots/information are desired. 

ques5=input('Would you like to see the final estimate of the Inverse Moment of 
Inertia Matrix?? If the answer is yes then type 0, otherwise type 1 .»>’); 
if ques5==0, 

Ii 

pause; 

Eg 

pause; 

clg; 

subplot(221); 
plot(t,est_x(l ,:)),grid 

title(' Inverse Moment of Inertia Estimate- 111’) 

xlabel('Time - seconds') 

subplot(222); 

plot(t,est_x(2,:)),grid,title(’I 1 2') 
xlabel(’Time - seconds’) 
subplot(223); 

plot(t,est_x(3,:)),grid,title(’1 1 3') 
xlabel('Time - seconds’) 
subplot(224); 

plot(t,est_x(4,:)),grid,title(’I22’) 
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xlabel('Time - seconds'), pause; 
clg; 

subplot(221); 

plot(t,est_x(5,:)),grid,xlabel(’Time - seconds') 
title(’Inverse Moment of Inertia Estimate - 123’) 
subplot(222); 

plot(t,est_x(6,:)),grid,xlabel('Time - seconds') 

title('I33’), pause; 

clg; 

plot(eigcheck(l,:)),grid,title('First Eigenvalue of the Estimated Inverse MOI 
Tensor’), pause 

plot(eigcheck(2, :)),grid,title('Second Eigen value'), pause 
plot(eigcheck(3,:)),grid,title('Third Eigenvalue’),pause 
% 

else, 

end 

% End of Program 



% PROGRAM NAME: CERSlewKl 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: guessmoi, actualmoi, EULERa2p, 

% EULERp2a321 

% 

% DESCRIPTION: 

% The below program will simulate the CER Control System for Large 
% Angle (Slewing) Motion. The Control Law is a Quaternion Feedback 
% Regulator. However, it is an Adaptive Controller in that the MOI 
% of the CER + Target is unknown and potentially variable. A KALMAN 
% Filter is used to estimate the MOI and then update the feedback gains 

% produced by the Quaternion Feedback Regulator. This KALMAN Filter 

% assumes that angular acceleration measurements are available, 
clear 

% Prompt user for first guess inverse MOI tensor 
% 

guessmoi; 

% Define the initial state estimate for the Kalman Filter 
est_x=[Ig(l,l);Ig(l,2);Ig(l,3);Ig(2,2);Ig(2,3);Ig(3,3)]; 

% 

% Define the actual MOI tensor. 

% This is placed into the MATLAB workspace by a separate subroutine, 
actualmoi; 

% 

T=75e-3; % The sampling time. 
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t(l)=l; %Initial sampling index 
% Initial orientation 
anginit=[0 0 0]; 



% 

xEULER(l :3, 1 )=[anginit( 1 );anginit(2);anginit(3)] ; 

% Prompt the user for the desired orientation 

angc=input('Specify the final or desired orientation in terms of 3-2-1 Euler 
Angles, in degrees. Enter as follows, [120 120 120]>»»'); 

% 

% Now, convert each of these orientations to Euler Parameters 
% Initial Conditions 

binit=EULERa2p(2,0,anginit(l),anginit(2),anginit(3)); 

beta(:,l)=[binit(l);binit(2);binit(3);binit(4)]; 

w(:,l)=[0;0;0]; 

% 

bcmd= EULERa2p(2,0,angc( 1 ),angc(2),angc(3)); 

% Compute the commanded quaternion matrix 

b0c=bcmd(l); 

blc=bcmd(2); 

b2c=bcmd(3); 

b3c=bcmd(4); 

BCMD=[b0c blc b2c b3c;-blc bOc b3c -b2c;-b2c -b3c bOc blc;-b3c b2c -blc 
bOc]; 

% 

% Prompt the user for the values used to compute the K and D weighting 
% matrices used in the Quaternion Feedback Regulator. 
settle=input(’Please enter the desired Settling Time in seconds»»'); 

% 

runtime=input(’Enter desired simulation run time (sec)»»'); 
actual=runtime/T; 

NUM=actual+0.2*actual ; 

NUM=round(NUM) ; 

% 

% The two weighting matrices, K and D will now be calculated as follows: 
zeta=1.0 

% 

omegaN=8/(zeta*settle); 
dscalar=2*zeta*omegaN ; 
kscalar=2*omegaN A 2; 

% 

rand('normar); 

% 

rd=(1.4e-4)*rand(3,NUM); % Random Plant Disturbance 



69 



mnl=((1.745e-4)/3)*rand(4,NUM); % Random Measurement Noise 
mn2=((5e-6)/3)*rand(3,NUM);% 

% 

% Initial Kalman Filter Parameters 
phik=eye(6); % Define the phi matrix for the plant 
% Prompt user for Kalman Filter Design Plant Noise 

Wnoise=input('Specify the scalar multiplier for the Kalman Filter Plant Noise 
Matrix. < 0 > Default Value (TRUST ME). < 1 > No NOISE. < 2 > Choose 
your own value.>>>>>'); 
if Wnoise==0, 

W=le-8; % The scalar multiplier of the below matrix 

% 

elseif Wnoise==l, 

W=0; 

else, 

W=input(’Specify your own value.»>»’); 
end 

Ew=W*eye(6); % Plant Noise Covariance Matrix 
V=1.0e-04; 

Ev=V*eye(3); % Measurement Noise Covariance Matrix 
% Prompt user again, for Covariance Estimation Error 
ques3=input('Specify the Covariance Estimation Error, P(0/0). < 0 > The 

Default. < 1 > The Identity Matrix. < 2 > A matrix of zeros.»>'); 

% 

if ques3==0, 

oldP=inv([-0.0114 -1.14e-6 2.92e-5 2.166e-5 4.84e-7 2.95e-5 

-1.14e-6 0.0001 -2.56e-7 -1.9e-7 -4.25e-9 -2.587e-7 

2.92e-5 -2.56e-7 -0.0026 4.864e-6 1.087e-7 6.623e-6 

le-6 le-6 le-6 -.0019 le-6 le-6 

le-6 le-6 le-6 le-6 -4.246e-5 le-6 

le-6 le-6 le-6 le-6 le-6 -2.587e-3]); 

elseif ques3==l, 

oldP=eye(6); 

else, 

oldP=zeros(6,6); 

end 

% This is P(0/0), the covariance matrix of the estimation error 

% 

tclock0=clock; 

% Initialize omega with noise 
wN(:,l)=w(:,l); 

% Now, simulate the CER + Target 
for N=1:NUM 
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% 

wN(:,N+l)=w(:,N)+mn2(:,N); % angular velocity plus measurement noise. 

% 

betaN(:,N)=beta(:,N)+mnl(:,N); % Position measurements plus noise. 
wskewN=[0 -wN(3,N) wN(2,N);wN(3,N) 0 -wN(l,N);-wN(2,N) wN(l,N) 0]; 

% 

wskew=[0 -w(3,N) w(2,N);w(3,N) 0 -w(l,N);-w(2,N) w(l,N) 0]; 

% 

Gskew=.5*[0 -w(l,N) -w(2,N) -w(3,N);w(l,N) 0 w(3,N) -w(2,N);w(2,N) 
-w(3,N) 0 w(l,N);w(3,N) w(2,N) -w(l,N) 0]; 

% 

t(N+l)=N+l; % Increment time index 

% 

% 

% Define the error quaternion for the feedback equation 

be=BCMD*betaN(:,N); 

qe=[be(2);be(3);be(4)] ; 

% Define the weighting matrices K and D. Note that they are updated by the 
% Kalman Filter’s Estimate of the MOI. 

K=kscalar*Ig; 

Dq=dscalar*Ig; 

U(:,N)=wskewN*Ig*wN(:,N) - Dq*wN(:,N) - K*qe; 

% Now, apply the limits for the control torques 

if abs(U(l,N» > 4 

U( 1 ,N)=4* sign(U( 1 ,N)); 

else 

U(1,N)=U(1,N); 

end 

if abs(U(2,N)) > 3 

U(2,N)=3*sign(U(2,N)); 

else 

U(2,N)=U(2,N); 

end 

if abs(U(3,N)) > 3 

U(3,N)=3*sign(U(3,N)); 

else 

U(3,N)=U(3,N); 

end 

% 

torque(:,N)=U(:,N)+rd(:,N); % Apply the random disturbance torques 
% 

% Discrete Nonlinear Kinematic and Dynamic Equations 

w(:,N+l)=w(:,N)+T*(Ii*torque(:,N)-Ii*wskew*MOI*w(:,N)); 
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beta(:,N+l)=beta(:,N)+T*Gskew*beta(:,N); 

% The real angular acceleration 
wdotreal(:,N)=Ii*torque(:,N)-Ii*wskew*MOI*w(:,N); 

% Angular acceleration from numerical differentiation 
wdot(:,N)=(wN(:,N+l)-wN(:,N))./T; 

% 

% Now, call the Kalman Filter in order to estimate the 
% actual MOI and update the Control Law. 

% 



Ck=[wdot( 1 ,N) (wdot(2,N)-wN(l,N)*wN(3,N)) (wN(l ,N)* wN(2,N)+wdot(3,N)) 
(-wN(2,N)*wN(3,N)) (wN(2,N) A 2-wN(3,N) A 2) (wN(2,N)*wN(3,N)) 
wN(l,N)*wN(3,N) (wdot(l,N)+wN(2,N)*wN(3,N)) (wN(3,N) A 2-wN(l,N) A 2) 
wdot(2,N) (wdot(3,N)-wN( 1 ,N)*wN(2,N)) >wN(l,N)*wN(3,N) 

-wN(2,N)*wN(l ,N) (wN(l ,N) A 2-wN(2,N) A 2) (wdot(l ,N)-wN(2,N)*wN(3,N)) 
wN(l,N)*wN(2,N) (wdot(2,N)+wN(l,N)*wN(3,N)) wdot(3,N)]; 

% 

% Below are the Kalman Filter and gain equations 

% 

% First, the Gain Equations 

newP=phik*oldP*phik' + Ew; 

G=newP*Ck , *inv(Ck*newP*Ck’ + Ev); 
oldP=(eye(6) - G*Ck)*newP; 

% 

est_xx=phik*est_x(:,N); 
est_y= Ck*est_xx; 

est_x(:,N+l)= est_xx + G*(torque(:,N)-est_y); 

% 

II l=est_x(l,N); 

I12=est_x(2,N); 

I13=est_x(3,N); 

I22=est_x(4,N); 

I23=est_x(5,N); 

I33=est_x(6,N); 

IEST=[I1 1 112 I13;I12 122 I23;I13 123 133]; 

% 

% Check the inverse MOI guess for physical reality. Eigenvalues must be > 0. 
eigcheck(l :3,N)=[eig(EEST)]; 

% 

checke(N)=l *sign(eigcheck(l ,N))+1 *sign(eigcheck(2,N))+l *sign(eigcheck(3,N)) 



% 

if checke(N)==3, 
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Ig=IEST; 



Update the initial guess MOI 
Otherwise, do not update the control law 



else, 

end 

% 

% 

% Convert the Euler Parameters to Euler Angles 

xEULER( : , N+ 1 )=EULER P 2a3 2 1 (beta( : , N+ 1 )) ; 
end % End of Simulation Loop 

% 

% Convert the time vector to actual time 
t=t.*T; 

% create a special time vector for the control input plot 

tU=0:N-l; 

tU=tU.*T; 

etime(clock,tclock0)/60 

% Calculate the total fuel used (actually this is proportional to it) 

FUELUSED=sum(sum(abs(U))) 

% 

% End of CER Slewing Simulation Program 



% PROGRAM NAME: replot 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: plot, subplot 

% 

% DESCRIPTION: 

% Plotting program for previously calculated data 

dgj 

subplot(221); 

plot(t,xEULER(l,:)),grid 

xlabel('Time - seconds'), ylabel(’ Angle - degrees') 

title('Z or Yaw Angle History’) 

% 

subplot(222); 
plot(t,xEULER(2, :)),grid 

xlabel(’Time - seconds'), ylabel(' Angle - degrees') 
title('Y or Pitch Angle History') 

% 

subplot(223); 

plot(t,xEULER(3,:)),grid 

xlabel('Time - seconds'), ylabel('Angle - degrees') 

title('X or Roll Angle History’) 
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% 

subplot(224); 

plot(t,w(l,:),t,w(2,:),’--’,t,w(3, :),’-.’), grid 

xlabel('Time - seconds’), ylabel('Angular velocity - rad/sec') 

title('Angular Velocity History’) 

pause; 

clg; 

% 

subplot(221); 

plot(beta(2,:),beta(3,:)),grid 

title(’Quatemion Plot- q2 vs. qr),xlabel('ql’),ylabel(’q2’) 

% 

subplot(222); 

plot(beta(2,:),beta(4,:)),grid 

title(’Quatemion Plot- q3 vs. qr),xlabel('qr),ylabel(’q3’) 

% 

subplot(223); 

plot(beta(3,:),beta(4,:)),grid 

title('Quatemion Plot- q3 vs. q2’),xlabel('q2’),ylabel(’q3’) 

% 

subplot(224); 

plot(tU,U(l,:),tU,U(2, tU,U(3, :),’-.’), grid 
title(’Control Torque History') 

xlabel('Time - seconds’), ylabel('Control Input/Torque - foot-lbs’) 

pause; 

clg; 



% Prompt user to see if additional plots/information are desired. 
ques5=input('Would you like to see the final estimate of the Inverse Moment of 
Inertia Matrix?? < 0 > YES. < 1 > NO »>'); 
if ques5==0, 

MOI 

pause; 

MOIestimate=IEST 

pause; 

invMOI=Ii 

pause; 

invMOIestimate=inv(IEST) 

pause; 

MOIeigenvalues=eig(MOI) 

pause; 

MOIesteigenvalues=eig(MOIestimate) 

pause; 
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subplot(221); 

plot(t,est_x(l,:)),grid 

title('Moment of Inertia Estimate - 111') 

xlabel('Time - seconds’) 

subplot(222); 

plot(t,est_x(2,:)),grid,title('I12’) 
xlabel(’Time - seconds’) 
subplot(223); 

plot(t,est_x(3,:)),grid,title('Il 3’) 
xlabel(’Time - seconds’) 
subplot(224); 

plot(t,est_x(4,:)),grid,title(’I22’) 
xlabel(’Time - seconds'), pause; 
clg; 

subplot(221); 

plot(t,est_x(5,:)),grid,xlabel('Time - seconds’) 
title('Moment of Inertia Estimate - 123’) 
subplot(222); 

plot(t,est_x(6,:)),grid,xlabel(’Time - seconds') 

title(’I33’) 

subplot(223); 

% Plot fuel usage plot 
plot(sum(abs(U))),grid,title(’Fuel Usage’), pause; 

% 

clg; 

% 

plot(eigcheck(l ,:)),grid,title(’First Eigenvalue of the Estimated MOI 
Tensor’), pause 

plot(eigcheck(2,:)),grid,title(’Second Eigenvalue'), pause 
plot(eigcheck(3,:)),grid,title(’Third Eigenvalue'), pause 
% 

else, 

end 

% 

% End of Program 



% PROGRAM NAME: EULERa2p 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% The below function will compute the four Euler parameters 

% given a sequence of three Euler Angles. A 3-2-1 or a 3-1-2 
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% rotation sequence is allowed. 

% 

% function E=EULERa2p(e,r,al,a2,a3); 

% 

% Here are the required input arguments: 

% Is the sequence 3-2-1 (e=2) or 3-1-2 (e=l)? 

% Are the three angles in radians ? (Yes: r=l, No: r=0) 

% The three Euler Angles (al,a2,a3) 

% 

% A vector is returned that contains b0,bl,b2,b3. 
function E=EULERa2p(e,r,al ,a2,a3); 

% 

% Convert the angles to radians if necessary 

if r==0 

% 

al=al*pi/180; 

a2=a2*pi/180; 

a3=a3*pi/180; 

% 

else 

end 

% Define the R matrices 

R0=[1 0 0 0;0 1 0 0;0 0 1 0;0 0 0 1]; 

R1=[0 -1 0 0;1 0 0 0;0 0 0 1;0 0 -1 0]; 

R2=[0 0 -1 0;0 0 0 -1;1 0 0 0;0 1 0 0]; 

R3=[0 0 0 -1;0 0 1 0;0 -1 0 0;1 0 0 0]; 

% 

% Test to see which rotation sequence is desired and proceed. 

% 

if e==l % This is a 3-1-2 sequence 
R3 1 2=(cos(a3/2)*R0 

+sin(a3/2)*R2)*(cos(a2/2)*R0+sin(a2/2)*Rl)*(cos(al/2)*R0+sin(al/2)*R3); 

E=R312*[1;0;0;0]; 

else 

% This is a 3-2-1 sequence. 

R321=(cos(a3/2)*R0+sin(a3/2)*Rl)*(cos(a2/2)*R0+sin(a2/2)*R2)*(cos(al/2)*R0 

+sin(al/2)*R3); 

E=R321*[1;0;0;0]; 

end 

% 

% End of function 



[% PROGRAM NAME: EULERp2a321 
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% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% The below function will calculate a sequence of three Euler Angles 
% from given Euler Parameters. A 3-2-1 rotation sequence is 
% allowed. 

% 

% function E=EULERp2a321(bvector) 

% 

% Here is the required input argument: 

% The 4 Euler Parameters: b0,bl,b2,b3 

% Enter these parameters as a vector. For example, 

% bvector=[b0;bl;b2;b3] 

% A vector is returned that contains the three angles in degrees, 
function E=EULERp2a321(bvector); 

% 

bO=bvector(l); 

bl=bvector(2); 

b2=bvector(3); 

b3=bvector(4); 

% First, compute the Direction Cosine Matrix (DCM) from the Euler 
Parameters. 

Cp=[bO A 2+bl A 2-b2 A 2-b3 A 2 2*(bl *b2+b0*b3) 2*(bl*b3-b0*b2);2*(bl*b2- 
b0*b3) bO A 2-b 1 A 2+b2 A 2-b3 A 2 2*(b2*b3+b0*bl);2*(bl*b3+b0*b2) 2*(b2*b3- 
bO*bl) b0 A 2-bl A 2-b2 A 2+b3 A 2]; 

% 

% Calculate the middle angle, beta 

a2=atan2(-Cp( 1 ,3 ),((Cp( 1 , 1 )) A 2+(Cp( 1 ,2)) A 2) A .5); 

% Test this angle and calculate the other angles 
a2d=a2*( 180/pi); 
a2d=round(a2d); 
if a2d==90, 

al=0;% This is the angle alpha or the rotation about the z-axis 
a3=atan2(Cp(2, 1 ),Cp(2,2)); 

% 

elseif a2d==-90, 
a 1=0; 

a3=-atan2(Cp(2, 1 ),Cp(2,2)); 

% 

else, 

al=atan2((Cp( 1 ,2)/cos(a2)),(Cp( 1 , 1 )/cos(a2»); 
a3=atan2((Cp(2,3)/cos(a2)),(Cp(3,3)/cos(a2))); 
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end 



% 

E=[al;a2;a3]; 

E=(180/pi)*E; 

% 

% End of function 



% PROGRAM NAME: Inew 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% The below function will calculate the new Moment of Inertia (MOI) 
% Tensor for the CER plus the Target. 

% The Target's MOI can be in the CER ref frame or its own; 

% as long as a Direction Cosine Matrix (DCM) is available 
% function I=Inew(Icer,Itar,Ml,M2,r2,ref,dcm); 

% 

% Here are the required input arguments: 

% MOI of CER (leer) 

% MOI of Target (Itar) 

% ref (If it equals 0 then Itar is in the CER ref frame) 

% DCM from the Target ref frame to the CER ref frame (dem) 

% Mass of the CER (Ml) 

% Mass of the Target (M2) 

% The vector between the Center of Mass of the CER 
% and the target mass in Cartesian Coord. (r2) 

% 

function I=Inew(Icer,Itar,Ml ,M2,r2,ref,dcm); 

% First, test to see if the Target MOI is in the CER ref frame 
if ref==0 
Itar = Itar; 
else 

Itar= dcm*Itar*dcm'; 
end 

% 

% Define the elements of the R2 Matrix 
R2=[0 -r2(3) r2(2);r2(3) 0 -r2(l);-r2(2) r2(l) 0]; 

MF=M1*M2/(M1+M2); % The Mass Factor multiplying R2 

% 

1= leer + Itar - MF*R2*R2; 

% End of Function 
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% PROGRAM NAME: DCM 

% PROGRAM AUTHOR: Nicholas F. Russo 

% CALLED PROGRAMS: none 

% 

% DESCRIPTION: 

% The below function will calculate the Direction Cosine Matrix 
% (DCM) for a 3 Euler Angle rotation sequence. A 3-2-1 (Yaw, 
% Pitch, Roll) or a 3-1-2 (Yaw, Roll, Pitch) is allowed. 

% 

% function D=DCM(e,r,al,a2,a3); 

% 

% Here are the required input arguments: 

% Is the sequence 3-2-1 (e=2) or 3-1-2 (e=l) ? 

% The three Euler Angles (al,a2,a3), in radians. 

% Are the Euler Angles in radians (Yes: r=l, No: r=0) 

function D=DCM(e,r,al,a2,a3); 

% 

% Convert the angles to radians if necessary 

if r==0 

al=al*pi/180; 

a2=a2*pi/l 80; 

a3=a3*pi/180; 

else 

end 

% First, set up the individual rotation matrices. 

C3=[cos(al) sin(al) 0;-sin(al) cos(al) 0;0 0 1]; 

% This corresponds to a rotation about the # 3 axis 
% 

C2a=[cos(a2) 0 -sin(a2);0 1 0;sin(a2) 0 cos(a2)] ;% # 2 axis 
% 

C2b=[cos(a3) 0 -sin(a3);0 1 0;sin(a3) 0 cos(a3)]; % # 2 axis 
% 

Cla=[l 0 0;0 cos(a3) sin(a3);0 -sin(a3) cos(a3)]; % # 1 axis 
% 

Clb=[l 0 0;0 cos(a2) sin(a2);0 -sin(a2) cos(a2)]; % # 1 axis 
% 

% 

% Test to see which rotation sequence is desired and proceed. 

% 

if e==l 

D= C2b*Clb*C3; 
else 

D= Cla*C2a*C3; 
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D= Cla*C2a*C3; 
end 

% 

% End of this function 
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